import numpy as np
import matplotlib.pyplot as plt
from math import pi
def chebypolys1st(t,order):
    deg=order+1
    M=np.zeros((t.size,deg))
    M[:,0]=1
    M[:,1]=t.T
    for n in range(2,deg):
        M[:,n]=2*t*M[:,n-1]-M[:,n-2]
    return M

def chebyfit1st(t,x,order,rcond=1e-10,coef=True):
    M=chebypolys1st(t,order)
    c=np.linalg.lstsq(M,x.T,rcond)[0]
    c=c.reshape(1,c.size)
    if coef:
        return c.reshape(c.size)
    else:
        f=np.dot(M,c.T)
        return f

def cheby1stval(t,c):
    M=chebypolys1st(t,c.size-1)
    c=c.reshape(1,c.size)
    f=np.dot(M,c.T)
    return f

def hilbertcheby(x,order,taper=0,rcond=1e-10,coef=True):
    L=x.size
    d=2.0/(x.size-1)
    t=np.arange(-1,1+d,d)
    if taper!=0:
	if taper>1:
		print "taper should be less than 1!"
		return -1
	else:
		tpz=int(L*taper/2.0)
		scl=np.ones(x.size)
		scl[0:tpz]=(1-np.cos(pi/tpz*np.arange(0,tpz)))*0.5
		scl[x.size-tpz:]=(1+np.cos(pi/tpz*np.arange(0,tpz)))*0.5
		sx=scl*x
    else:
	sx=x			
    a=chebyfit1st(t,sx,order,rcond=rcond)			
    c=np.zeros(a.size)
    for i in range(0,c.size):
        for j in range(0,a.size):
            if i>j:
                c[i]=c[i]+1/pi*a[j]*((1-np.power(-1,i+j+1))/(i+j+1.0)+(1-np.power(-1,i-j+1))/(i-j+1.0))
            else:
                c[i]=c[i]+1/pi*a[j]*(1-np.power(-1,i+j+1))/(i+j+1.0)
    if coef:
        return c
    else:
        hf=cheby1stval(t,c)
        return hf

def dhilbertcheby(x,order,taper=0,rcond=1e-10):
    L=x.size
    d=2.0/(x.size-1)
    t=np.arange(-1,1+d,d)
    c=hilbertcheby(x,order,taper,rcond)
    U=np.zeros((L,c.size))
    U[:,0]=1
    U[:,1]=2*t.T
    for i in range(2,c.size):
        U[:,i]=2*t*U[:,i-1]-U[:,i-2]
    s=np.arange(0,c.size)+np.ones(c.size)
    sc=s*c
    f=np.dot(U,sc.T)
    return f

if __name__=="__main__":
    t=np.arange(0,256)
    x=np.sin(pi/20.0*t)
    y=dhilbertcheby(x,50,taper=0.05)
    plt.plot(x)
    plt.plot(y,'r')
    plt.show()
    
